function [g_hat] = npreg(X,Y)
%NPREG nonparametric regression using kernel

[n,q] = size(X);

h_vec = std(X)*n^(-1/(4+q));  % rule-of-thumb bandwidth
kernel = @(x)normpdf(x); % normal kernel

g_hat = zeros(n,1);
for i = 1 : n
    K_i = prod(kernel((X-repmat(X(i,:),n,1))./repmat(h_vec,n,1)),2);
    g_hat(i) = dot(Y,K_i)/sum(K_i);
end

